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We study the fragility of spin glasses to small temperature perturbations numerically using population an¬ 
nealing Monte Carlo. We apply thermal boundary conditions to a three-dimensional Edwards-Anderson Ising 
spin glass. In thermal boundary conditions all eight combinations of periodic versus antiperiodic boundary con¬ 
ditions in the three spatial directions are present, each appearing in the ensemble with its respective statistical 
weight determined by its free energy. We show that temperature chaos is revealed in the statistics of crossings 
in the free energy for different boundary conditions. By studying the energy difference between boundary con¬ 
ditions at free-energy crossings, we determine the domain-wall fractal dimension. Similarly, by studying the 
number of crossings, we determine the chaos exponent. Our results also show that computational hardness in 
spin glasses and the presence of chaos are closely related. 


Chaos refers to sensitivity to small perturbations. In ad¬ 
dition to dynamical systems where the phenomenon was 
first identified, there are many statistical mechanical systems 
where chaotic effects have been predicted and observed. For 
example, hysteresis, memory, and rejuvenation effects found 
in random elastic manifolds, polymers CM, as well as spin 
glasses are considered to be a direct manifestation of the pres¬ 
ence of chaos isa. It is surprising and fascinating that both 
the nonequilibrium and equilibrium states of spin glasses are 
so fragile to small perturbations. Chaos is therefore central 
to the understanding of both equilibrium and nonequilibrium 
properties of spin glasses, as well as related systems. The con¬ 
nection between chaos in spin glasses and dynamical systems 
has been recently explored [|8l. Furthermore, there is mount¬ 
ing evidence that chaos in spin glasses is directly related to 
the computational hardness and long thermalization times 0 
of these paradigmatic benchmark problems. As such, quan¬ 
tifying and understanding chaotic effects in spin-glass-like 
Hamiltonians could be of great importance for the develop¬ 
ment of any novel algorithm or computing architecture Ha- 

ini. 

In this work we study the effects of thermal perturbations. 
Temperature chaos refers to the property that a small change 
in temperature results in a complete reorganization of the 
equilibrium configuration of the system. Temperature chaos 
has long been predicted for spin glasses lfT3l - [T^ . Although 
some early studies raised doubts about the existence of tem¬ 
perature chaos jm, increasing numerical evidence for tem¬ 
perature chaos has emerged in recent years for various mod¬ 
els such as the random-energy random-entropy model llTSi 
and also more realistic three- and four-dimensional Ising spin 
glasses 011211201. It has been suggested that temperature 
chaos would only be observable in spin glasses at very large 
system sizes and large changes in the temperature 11211 12^ . 
However, some studies EOl demonstrated the existence of 
temperature chaos via scaling arguments. 

One direct manifestation of temperature chaos is that the 
free-energy difference between two boundary conditions that 
differ by a domain wall may change sign as a function of tem¬ 


perature. Previous studies examined the free-energy differ¬ 
ence between periodic and antiperiodic boundary conditions 
in a single direction to identify temperature chaos nuEii. 
This motivates us to study temperature chaos using thermal 
boundary conditions in which all 2“^ combinations of pe¬ 
riodic and anti-periodic boundary conditions in the d spatial 
directions appear in a single simulation with their appropri¬ 
ate statistical weights. Thermal boundary conditions provide 
a novel and elegant way to study temperature chaos. 

Here we quantitatively investigate temperature chaos using 
population annealing Monte Carlo Il24lj28l . This simulation 
approach is ideal to study chaos effects in spin glasses because 
multiple boundary conditions can be studied at the same time. 
We show that temperature chaos is revealed in the statistics 
of crossings in the free energy for pairs of boundary condi¬ 
tions ll23l and thus establish both qualitatively and quantita¬ 
tively the presence of chaos in spin glasses. Our approach can 
be applied to a multitude of problems and, in particular, to 
the search for hard benchmark instances for novel computing 
paradigms ifTTlflEl . 

What causes temperature chaos? Temperature chaos re¬ 
sults from the existence of dissimilar classes of configura¬ 
tions with similar free energies but differing energies and 
entropies Gsiiia. Consider two classes of spin configura¬ 
tions, a I and CT 2 , corresponding to distinct basins in the free- 
energy landscape. Within each class, all spin configurations 
are similar but the two classes are dissimilar and differ by a 
large relative domain wall. Let AF{T) be the free-energy 
difference at temperature T between these two classes, with 
AF{T) = AE{T) - TAS{T) where AE and AS' are the 
energy and entropy, respectively, of the relative domain wall. 
Suppose now that AE and AS are both much larger than AE 
and weakly dependent on temperature; then a small change 
in temperature may lead to sign change in AE. Suppose 
that AF, AE, and AS all behave as power laws in the size 
scale i of the relative domain wall separating spin configu¬ 
rations CTi and CT 2 with leading behavior AF ^ but with 
AE ~ AS ~ and ds/2 > 9. Here 9 is the stiffness 
exponent and dg is the fractal dimension of the domain wall. 
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As £ increases, the temperature perturbation 5T required to 
change the sign of A_F decreases, i.e., 5T ~ with the 
chaos exponent C, given hy (^ = ds/2 — 6 ifThl . 

We investigate temperature chaos in the Edwards-Anderson 
(EA) Ising spin-glass model 1291 . The EA Hamiltonian is 

Td — ^ ^ JijSiSj ^ ( 1 ) 

{iJ) 

where Si = ±1 are Ising spins. The sum (ij) is over the 
nearest-neighbor sites in a cubic lattice with N = sites. 
Jij is the interaction between spins Si and sj, and is chosen 
from a Gaussian distribution with mean zero and variance 1. 
We refer to each disorder realization as a “sample.” 

We use thermal boundary conditions (TBC) to study tem¬ 
perature chaos in the EA model. In the TBC ensemble each 
boundary condition i occurs in the ensemble with a weight de¬ 
pending on its free energy Fi. The probability pi of boundary 
condition z in the ensemble is given by Pi = exp[—/3(Fi—F)], 
where F is the total free energy of the system in TBC and /3 
the inverse temperature. Thermal boundary conditions were 
introduced to minimize the finite-size effects due to domain 
walls and have proved to be useful in studying the low- 
temperature phase of the EA model l24l . They have been 
used with exact algorithms for finding ground states of two- 
dimensional spin glasses lEolliIl (referred to there as “ex¬ 
tended” boundary conditions). A more restricted version of 
TBC using periodic and antiperiodic boundary conditions in 
only a single direction was used in Refs. umiiiHii. 

In thermal boundary conditions, a domain wall on the scale 
of the linear system size L separates each boundary condition. 
Thus temperature chaos manifests itself as a strong tempera¬ 
ture dependence in the relative free energies of the different 
boundary conditions (BCs). Because the stiffness exponent 
is positive, in the low-temperature phase one expects that for 
large systems a single BC will dominate the ensemble for al¬ 
most all temperatures. However, as the temperature changes, 
the dominant boundary condition will frequently change. A 
crossing event occurs when the free-energy difference be¬ 
tween two BCs changes sign. The proliferation of crossing 
events is a direct indication of temperature chaos. Boundary- 
condition crossing events between periodic and antiperiodic 
BCs in one direction were studied in the two-dimensional EA 
model in Ref. ||2^ and identified as a signature of tempera¬ 
ture chaos. Eigure [T] shows BC probabilities pi for all eight 
boundary conditions as a function of temperature for a single 
L = 10 sample. As expected, at high temperatures, each BC 
occurs with equal probability. However, at low temperatures, 
four different BCs dominate in different temperature ranges 
and, indeed, the dominant boundary condition at the lowest 
temperatures has a tiny probability in a range just below the 
critical temperature. 

We carried out simulations of the three-dimensional EA 
model in TBC using population annealing Monte Carlo 
1^ . Population annealing is similar to simulated annealing: 
In both algorithms, the system is cooled from a high temper¬ 
ature to a low temperature following an annealing schedule. 
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FIG. 1: (Color online) A single size L = 10 sample displaying sev¬ 
eral boundary-condition crossing. The plot shows the probabilities 
of the eight boundary conditions {pi } as a function of inverse tem¬ 
perature j3. 


However, population annealing involves cooling a population 
of replicas and includes a resampling of the population as it 
is cooled. At each temperature step in the annealing sched¬ 
ule, each replica is acted on independently by the Metropolis 
algorithm. In the resampling step, which occurs before the 
temperature is changed, replica i is differentially reproduced 
according to its energy Ei. The expected number of copies 
of replica i is exp[—(^' — f])Ei]/Q{/3, /3') for a temperature 
step from /3 to /3'. The normalization Q{/3, P') is chosen such 
that the expected population size is unchanged by the resam¬ 
pling step, Q{I3,13') = (l/i?o) exp[-(/3' - /3)F,], where 

i?o is the expected population size. The actual number of 
copies made of replica i is a random integer whose mean is 
exp[—(/3' — f3)Ei\/Q{/3, P'). Note that expected number of 
copies of each replica is exactly the reweighting factor be¬ 
tween Gibbs distributions at /3 and (3'. Thus, if the population 
is representative of the Gibbs distribution at inverse temper¬ 
ature (3 and i?o is large, then after resampling, the popula¬ 
tion is representative the Gibbs distribution at (3'. The anneal¬ 
ing schedule consists of Np temperature steps equally spaced 
in (3 with Ns = 10 Metropolis sweeps at each temperature. 
Thermal boundary conditions are easily simulated in popula¬ 
tion annealing by initializing the population at /3 = 0 with 
1/8 of the population in each of the eight BCs ||28l. Re¬ 
sampling takes care of making sure that at every temperature, 
each BC appears with the correct statistical weight. We study 
2000 samples of sizes L = 4 (Rq = 5 10^, Np = 101), 6 
(Rq = 2 10^, Nt = 101), 8 (i?o = 5 10^ Nt = 201), 10 
(Ro = 10®, Nt = 301), and 12 (Rq = 10®, Np = 301), 
down to temperature T = 0.33. The critical temperature is 
Tc ~ 0.951 llTSll so the simulations include temperatures that 
are deep within the low-temperature phase. Eor some hard 
samples ll^ we use larger population sizes. In the case of 
L = 12 approximately 300 samples needed to be run with up 
to a factor of 10 larger population sizes. 
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Using population annealing we carry out a quantitative 
study of boundary condition crossings. The temperature dif¬ 
ference between crossing scales as so that the num¬ 
ber of crossing Nc in a fixed temperature interval scales as 
Nc Also, at crossings, we have that Af = 0 so 

that ds/2 can be obtained from the scaling of the average of 
AE at crossings as function of L. Finally, in a previous study 
we measured the stiffness exponent in TBC. We defined the 
sample stiffness A as A = log[//(1 — /)] where / is the prob¬ 
ability of being in the dominant BC, i.e., / = maxi[{pi}]. We 
measured 9 as the scaling of the median of A with system size 
L. Thus, within TBC we can independently measure all three 
exponents 9, ds/2, and Q and verify the relation C, = dsl‘2, — 9. 

Crossings can be divided into two classes: Dominant cross¬ 
ings are those such that the two equal BC probabilities at the 
crossing are larger than all other BC probabilities. All other 
crossings are subdominant. For large systems, the BC prob¬ 
ability at a subdominant crossing is expected to be typically 
suppressed by a factor exp(L“®) relative to the dominant BC 
and thus be increasingly difficult to observe in TBC simula¬ 
tions. To avoid finite-size corrections in counting crossings, 
here we focus on dominant crossings. On the other hand, 
for measuring AE = TAS {AS the change in entropy) at 
crossings we do not expect a distinction between dominant 
and subdominant crossings and, to improve statistics, we use 
all crossing with > 0.05. 

Figure is a log-log plot (base 10) of Nc vs L where 
Nc counts dominant crossing in the range /3 G (1.5,3.0). 
A simple power-law fit Nc yields ^ = 0.96(5). All 

quoted error bars are one standard deviation statistical er¬ 
rors. To test the effect of temperature on this exponent, we 
also calculated () from two smaller temperature ranges. For 
P G (1.5, 2.0) we find C = 1.07(8), and from /3 G (2.0,3.0) 
we find C, = 0.85(9). For higher temperatures, critical fluc¬ 
tuations may contaminate the measurement of the chaos ex¬ 
ponent while for lower temperatures the number of crossings 
is suppressed by the smallness of the entropy. We note that 
there is a significant trend to a smaller value of ^ at lower 
temperatures. If one assumes that a single exponent holds 
throughout the low-temperature phase, this trend suggests sig¬ 
nificant temperature-dependent finite-size corrections. The in¬ 
set to Fig. 1^ shows a histogram of the number of crossings in 
the range /3 G (1.5, 3.0) with pi > 0.05 for size L = 12 as 
a function of inverse temperature and reveals that the num¬ 
ber of crossings decreases with temperature, consistent with 
the fact that the entropy decreases with temperature so that in¬ 
creasingly large temperature changes are required to change 
the free-energy difference between BCs. In the large-volume 
limit, the number of dominant crossings per sample is ex¬ 
pected to become infinite but for size L = 12 temperature 
chaos events are infrequent but not rare-there are on aver¬ 
age 0.86 crossings with pi > 0.05 per sample in the range 
/3 G (1.5, 3.0) and 0.33 dominant crossings per sample in 
the same temperature range. An advantage of using bound¬ 
ary condition crossings in TBC is that temperature chaos is 
not a rare event for accessible system sizes in contrast to over¬ 


lap correlations in a single boundary condition where chaotic 
effects are weak in most samples 0. 
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FIG. 2; (Color online) Number of dominant crossing in the range 
/3 G (1.5,3.0) vs size L, for L = 4, 6, 8, 10, and 12. The straight 
line is the best power law fit (see text). The inset is a histogram of the 
number of all crossings with pi > 0.05 with respect to /3 for system 
size L = 12. 

Figure|^is a log-log plot (base 10) of the median and mean 
of the absolute energy difference |Ai?| vs L at all crossings 
in the range f3 G (1.5,3.0) such that pi > 0.05. A sim¬ 
ple power-law fit for the mean yields |Ai?| ^ with 

ds/2 = 1.18(2) with the same result for the median. We 
again test the effect of the temperature range on ds/2 by di¬ 
viding the /? range into two intervals, f3 G (1.5, 2.0) and 
/3 G (2.0,3.0), from which we obtain the the results for the 
mean ds/2 = 1.14(2) and ds/2 — 1.26(3), respectively. 
There is a significant trend toward larger values at lower tem¬ 
peratures, suggesting temperature-dependent finite-size cor¬ 
rections. 

Our results for the three-dimensional EA model, ds/2 = 
1.17(2) and (/ = 0.96(5), are comparable but slightly smaller 
than previous work: For example, ds/2 = 1.29(1) was found 
in Ref. ll^ based on perturbations of the ground state, and 
ds/2 = 1.31(1) was found in Ref. OtI based on the variance 
of the link overlap. Recall that our result for ds/2 in the lower 
temperature range is 1.26(3), which is within error bars of 
these zero-temperature results and suggests large temperature- 
dependent finite-size corrections. Our result, (/ — 0.96(5), is 
somewhat smaller than (/ = 1.04 found in Ref. EOl from the 
spin overlap between different temperatures. Combined with 
our estimate of 0 = 0.27(2) ll24ll we find that the predicted re¬ 
lation C = ds/2 — 0 is reasonably well satisfied by our results. 

Temperature chaos partially explains why spin-glass simu¬ 
lations are computationally costly 0. All known efficient al¬ 
gorithms for equilibrating three-dimensional spin glasses rely 
on coupled simulations at many temperatures. Algorithms in 
this class include parallel tempering Monte Carlo ll^ . the 
Wang-Landau algorithm and the algorithm used in this 
work, population annealing ll25l |2^ . In these algorithms. 
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FIG. 3: (Color online) Log-log plot of the mean (red circles) and 
median (blue squares) energy difference between dominant boundary 
conditions at crossing in the range (3 € (1.5, 3.0) for L = 4, 6, 8, 
10, and 12. The straight line is the best power law fit (see text). Error 
bars are smaller than the symbols. 


fast mixing at high temperatures provides new configurations 
to the low-temperature simulations. Temperature chaos de¬ 
creases the effectiveness of these algorithms because the con¬ 
figurations supplied from higher temperatures are often rather 
different from the important configurations at lower temper¬ 
atures. In TBC, temperature chaos means that BCs that are 
important at high temperature are unimportant at low temper¬ 
ature. This phenomenon is evident in Fig.[2 One might worry 
that boundary conditions that should be important at low tem¬ 
perature are completely lost at higher temperatures so that the 
simulations do not reach the correct TBC equilibrium. To ver¬ 
ify that this is not the case, we performed an additional check 
of the equilibration of the TBC ensemble by re-doing several 
hundred of the hardest L = 12 samples using an order of mag¬ 
nitude larger population sizes in the simulation and we found 
no difference in the number of crossings for any sample. 

A direct measure of hardness for population annealing for 
a given sample is the characteristic family size p. In popula¬ 
tion annealing, most of the original population is eliminated 
by successive resampling steps and the final population is de¬ 
scended from a small subset of the initial population. Every 
member of the final population can be uniquely assigned to a 
“family” descended from some member of the initial popula¬ 
tion. Let rii be the fraction of the low-temperature population 
descended from replica i in the initial population, p is then 
defined as 


p = lim iio exp 

Rq—^OO 




( 2 ) 


where i?o is the population size in the simulation (n^ = 0 
there is no contribution to the sum). In practice, p is mea¬ 
sured using the large Rq of the simulation. Since there may 
be correlations between members of the same family, the pop¬ 
ulation size i?o in the simulation must be much larger than 


p to assure a large number of independent measurements and 
small statistical errors. Thus samples with the largest p require 
the most computational resources to simulate. We have shown 
ll24ll that p is also strongly correlated with the integrated auto¬ 
correlation time for parallel tempering Monte Carlo, measured 
in Ref. HOl so that the same conclusions are likely to hold for 
parallel tempering. Figure shows the disorder average of 
logj^o p vs L measured at /3 = 3 for two different classes of 
disorder samples. The Nc = 0 class has no temperature chaos 
events (crossings withp^ > 0.05) in the range G (1.5,3.0) 
while the JVc > 0 class has one or more temperature chaos 
events in the same range. The error bars are smaller than the 
data points and the plots show that p scales exponentially in 
L, p (X exp(L/f), for both classes but that the characteristic 
size scale i = 1.27(14) is significantly smaller for the chaotic 
samples than for nonchaotic samples where £ = 1.62(10). 
It is an interesting question whether temperature chaos slows 
down all algorithms for spin glasses, not just those that de¬ 
pend on coupling multiple temperatures. Recent studies have 
shown that muni computationally hard instances for clas¬ 
sical algorithms are also computationally hard for quantum 
annealing machines, like the D-Wave Two quantum annealer. 
As such, by measuring p for a given sample, we have a simple 
way to uniquely classify the complexity of a given instance. 
This means that our approach is of great importance in the de¬ 
velopment of hard problems to discern whether quantum an¬ 
nealing can outperform simulated annealing simulations (see, 
for example. Refs. ||4TH451 ). 



FIG. 4: (Color online) The average of the log (base 10) of the hard¬ 
ness p, measured at /3 = 3, vs size L for two classes of samples, 
those without crossing in the range j3 G (1.5, 3.0), Nc = 0 (red 
circles) and those with at least one crossing in that range, Nc > 0 
(blue squares). Error bars are smaller than the symbols. The insets 
are histograms of values of p for the two classes of samples. 


We have seen that thermal boundary conditions allow us to 
identify temperature chaos with boundary condition crossings 
and provide a tool for studying chaos quantitatively even for 
the small sizes accessible to simulations. It would be interest¬ 
ing to apply these ideas to other types of chaos in spin glasses 
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such as bond chaos. We have also established that temperature 
chaos is a significant determinant of computational hardness 
for multicanonical algorithms but it remains an open question 
as to whether temperature chaos is correlated with hardness 
for all algorithms. 
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